#include <assert.h>
#include "mcPublic.h"

#include "mcmcHistogram.h"
#define PI 3.1415926535897931

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
Histogram::Histogram()
{
	cut_off = false;

	pbins	= NULL;
	count	= 0;

	x0		= 0.0;
	d		= 0.0;
	mean	= 0.0;
	outof_lower_bound = 0.0;
	outof_upper_bound = 0.0;
}

Histogram::~Histogram()
{
	Free();
}

void Histogram::Free(void)
{
	if (pbins != NULL)
		delete []pbins;
	pbins = NULL;
	count	= 0;

	x0		= 0.0;
	d		= 0.0;
	mean	= 0.0;
	outof_lower_bound = 0.0;
	outof_upper_bound = 0.0;
}

bool Histogram::AllocateSpace(const int size)
{
	if (size!=count && size>0)
	{
		if (pbins != NULL)
			delete []pbins;
		pbins = new double[size];
		count = size;
	}

	return true;
}

// operations
Histogram& Histogram::operator=(const Histogram &a)
{
	if (&a != this)
	{
		AllocateSpace(a.count);

		if (a.count > 0)
			for (int i=0; i<count; i++)
				pbins[i] = a.pbins[i];

		x0		= a.x0;
		d		= a.d;
		mean	= a.mean;
		cut_off = a.cut_off;
		outof_lower_bound = a.outof_lower_bound;
		outof_upper_bound = a.outof_upper_bound;
	}
	return *this;
}

int Histogram::TotalCodingBits(void)
{
	int bits;

	bits = sizeof(float)+		//x0
		   sizeof(float);		//d
	for (int i=0; i<count; i++)
		bits += 2*sizeof(char);

	return bits;
}

Histogram& Histogram::operator+(Histogram &a)
{
	assert(count==a.count);
	assert(d==a.d);

	for (int i=0; i<count; i++)
		pbins[i] += a.pbins[i];

	outof_lower_bound += a.outof_lower_bound;
	outof_upper_bound += a.outof_upper_bound;

	return (*this);
}

double Histogram::operator*(Histogram &a)
{
	double	result;

	assert(count==a.count);
	assert(d==a.d);

	result = 0.0;
	for (int i=0; i<count; i++)
		result += pbins[i]*a.pbins[i];

	result += outof_lower_bound*a.outof_lower_bound;
	result += outof_upper_bound*a.outof_upper_bound;

	return result;
}

FILE & Histogram::operator>>(FILE &strm)
{
	fwrite(&count, sizeof(count), 1, &strm);
	fwrite(&x0, sizeof(x0), 1, &strm);
	fwrite(&d, sizeof(d), 1, &strm);
	fwrite(&mean,sizeof(mean),1,&strm);

	fwrite(&outof_lower_bound,sizeof(outof_lower_bound),1,&strm);
	fwrite(&outof_upper_bound,sizeof(outof_upper_bound),1,&strm);

	if (count>0)
		fwrite(&(pbins[0]),sizeof(pbins[0]), count, &strm);

	return strm;
}

FILE & Histogram::operator<<(FILE &strm)
{
	int count_temp;

	fread(&count_temp,sizeof(count_temp),1,&strm);
	AllocateSpace(count_temp);

	fread(&x0, sizeof(x0), 1, &strm);
	fread(&d, sizeof(d), 1, &strm);
	fread(&mean,sizeof(mean),1,&strm);

	fread(&outof_lower_bound,sizeof(outof_lower_bound),1,&strm);
	fread(&outof_upper_bound,sizeof(outof_upper_bound),1,&strm);

	if (count_temp>0)
		fread(&(pbins[0]),sizeof(pbins[0]), count_temp, &strm);

	return strm;
}

void Histogram::WriteToTxtFile(FILE &strm)
{
	fprintf(&strm,"%d\n",count);
	fprintf(&strm,"%.14lf\n",x0);
	fprintf(&strm,"%.14lf\n",d);
	fprintf(&strm,"%.14lf\n",outof_lower_bound);
	fprintf(&strm,"%.14lf\n",outof_upper_bound);

	for (int i=0; i<count; i++)
		fprintf(&strm, "%.14lf\n", pbins[i]);
}

void Histogram::ReadFromTxtFile(FILE &strm)
{
	int	count_temp;

	fscanf(&strm,"%d",&count_temp);
	fscanf(&strm,"%lf",&x0);
	fscanf(&strm,"%lf",&d);
	fscanf(&strm,"%lf",&outof_lower_bound);
	fscanf(&strm,"%lf",&outof_upper_bound);

	AllocateSpace(count_temp);
	for (int i=0; i<count_temp; i++)
		fscanf(&strm, "%lf", &(pbins[i]));
}


void Histogram::Set(const double x1, const double x2, const double d1)
{
	x0 = x1;
	d = d1;

	assert(x1<=x2);
	int count_temp = (int)floor((x2-x1)/d1)+1;
	assert(count_temp>0);

	AllocateSpace(count_temp);
	Init(0.0);
}

void Histogram::Init(const double a)
{
	for (int i=0; i<count; i++)
		pbins[i] = a;
	outof_lower_bound = 0.0;
	outof_upper_bound = 0.0;
}

// distance measure between two histograms
double Histogram::DistTV(Histogram &a)
{
	double	dist=0.0;

	assert(count==a.count);

	for (int i=0; i<count; i++)
		dist += fabs(pbins[i]-a.pbins[i]);
	return dist;
}

double Histogram::DistKL(Histogram &a)
{
	double	dist=0.0,p1,q1;

	assert(count==a.count);

	for (int i=0; i<count; i++)
	{
		p1 = max(pbins[i], 1e-7);
		q1 = max(a.pbins[i], 1e-7);

		dist += p1*log(p1/q1);
	}
	return fabs(dist);
}


double Histogram::Entropy(void)
{
	double entropy=0.0;

	for (int i=0; i<count; i++)
	{
		if (pbins[i]>0.0)
			entropy += -pbins[i]*log(pbins[i]);
	}
	return entropy;
}

double Histogram::DistToGaussian(void)
{
	double			*p_gauss;
	int				i;
	double			sigma,d,sum;

	// initialization
	p_gauss = new double[count];
	sigma = 0.0;
	mean = 0.0;

	// set the Gaussian
	for (i=0; i<count; i++)
		mean += pbins[i]*((double)i);
	for (i=0; i<count; i++)
		sigma += pbins[i]*((mean-i)*(mean-i));
	sigma = sqrt(sigma);

	sum = 0.0;
	for (i=0; i<count; i++)
	{
		p_gauss[i] = 1.0/(sqrt(2.0*PI)*sigma)*exp(-0.5*(pow((mean-i)/sigma,2.0)));
		sum += p_gauss[i];
	}
	if (sum>0.0)
		for (i=0; i<count; i++)
			p_gauss[i] = p_gauss[i]/sum;

	// compute distance
	d = 0.0;
	for (i=0; i<count; i++)
		d += fabs(p_gauss[i]-pbins[i]);

	return d;
}

double Histogram::ComputeMean(void)
{
	mean = 0.0;
	// set the Gaussian
	for (int i=0; i<count; i++)
		mean += pbins[i]*((double)i);

	return mean;
}

void Histogram::Normalize(void)
{
	double sum=0.0;

	for (int i=0; i<count; i++)
		sum += pbins[i];
	if (!cut_off)
	{
		sum += outof_lower_bound;
		sum += outof_upper_bound;
	}

	if (sum>0.0)
	{
		for (int i=0; i<count; i++)
			pbins[i] = pbins[i]/sum;
		if (!cut_off)
		{
			outof_lower_bound = outof_lower_bound/sum;
			outof_upper_bound = outof_upper_bound/sum;
		}
	}
}

mcmcVector<double> Histogram::BinDifference(Histogram &a)
{
	assert(count==a.count);

	mcmcVector<double> vc_dist;
	vc_dist.SetRowVector(count);
	for (int i=0; i<count; i++)
		vc_dist(i) = pbins[i]-a.pbins[i];
	return vc_dist;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
// Lagrange multipliers
Gibbs_B& Gibbs_B::operator=(Gibbs_B& a)
{
	if (&a != this)
		*((Histogram*)this) = *((Histogram*) &a);
	return *this;
}

FILE & operator>>(Gibbs_B &g, FILE &strm)
{
	*((Histogram *)&g)>>strm;

	return strm;
}

FILE & operator<<(Gibbs_B &g, FILE &strm)
{
	*((Histogram *)&g)<<strm;

	return strm;
}

// operations
double	Gibbs_B::Projection(Histogram &hist)
{
	assert(count==hist.count);
	assert(d==hist.d);

	double sum=0.0;
	for (int i=0; i<count; i++)
		sum += pbins[i]*hist.pbins[i];
	return sum;
}

double	Gibbs_B::ComputeZ(void)
{
	double Z;

	Z = 0.0;
	for (int i=0; i<count; i++)
		Z += exp(-pbins[i]);

	mean = 0.0;
	for (int i=0; i<count; i++)
		mean += pbins[i]*((double)i);
	mean = (x0+mean*d);
	return Z;
}


#undef PI
